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Spatially resolved genetic data is increasingly used to reconstruct the migrational history of 
species. To assist such inference, we study, by means of simulations and analytical methods, the 
dynamics of neutral gene frequencies in a population undergoing a continual range expansion in 
one dimension. During such a colonization period, lineages can fix at the wave front by means of a 
"surfing" mechanism [Edmonds C.A., Lillie A.S. & Cavalli-Sforza L.L. (2004) Proc Natl Acad Sci 
USA 101: 975-979]. We quantify this phenomenon in terms of (i) the spatial distribution of lineages 
that reach fixation and, closely related, (ii) the continual loss of genetic diversity (heterozygosity) 
at the wave front, characterizing the approach to fixation. Our simulations show that an effective 
population size can be assigned to the wave that controls the (observable) gradient in heterozygosity 
left behind the colonization process. This effective population size is markedly higher in pushed 
waves than in pulled waves, and increases only sub-linearly with deme size. To explain these and 
other findings, we develop a versatile analytical approach, based on the physics of reaction-diffusion 
systems, that yields simple predictions for any deterministic population dynamics. 



Population expansions in space are common events in 
the evolutionary history of many species [H, 0,0, 0,0, 0,0] 
and have a profound effect on their genealogy. It is widely 
appreciated that any range expansion leads to a reduc- 
tion of genetic diversity ("Founder Effect") because the 
gene pool for the new habitat is provided only by a small 
number of individuals, which happen to arrive in the un- 
explored territory first. In many species, the genetic foot- 
prints of these pioneers are still recognizable today and 
provide information about the migrational history of the 
species. For instance, a frequently observed south-north 
gradient in genetic diversity ( "southern richness to north- 
ern purity" [8[ ) on the northern hemisphere is thought to 
reflect the range expansions induced by the glacial cycles. 
In the case of humans, the genetic diversity decreases 
essentially linearly with increasing geographic distance 
from Africa 0,0], which is indicative of the human mi- 
gration out of Africa. It is hoped 0] , that the observed 
patterns of neutral genetic diversity can be used to infer 
details of the corresponding colonization pathways. 

Such an inference requires an understanding of how a 
colonization process generates a gradient in genetic diver- 
sity, and which parameters chiefly control the magnitude 
of this gradient. Traditional models of population genet- 
ics [13], which mainly focus on populations of constant 
size and distribution, apply to periods before and after a 
range expansion has occurred, when the population is at 
demographic equilibrium. However, the spatio-temporal 
dynamics in the transition period, on which we focus in 
this article, is less amenable to the standard analytical 
tools of population genetics, and has been so far stud- 
ied mostly by means of simulations [ll|, [l2|, [H, 0, EH . 
An analytical understanding is available only for a lin- 
ear stepping stone model in which demes (lattice sites) 
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are colonized one after the other, following deterministic 
logistic growth [l6[ or instantaneously 17], in terms of 
recurrence relations. 

Recent computer studies suggest that the neutral ge- 
netic patterns created by a propagating population wave 
might be understood in terms of the mechanism of "gene 
surfing" 0,0]: As compared to individuals in the wake, 
the pioneers at the colonization front are much more suc- 
cessful in passing their genes on to future generations, not 
only because their reproduction is unhampered by lim- 
ited resources but also because their progeny start out 
from a good position to keep up with the wave front (by 
means of mere diffusion) . The offspring of pioneers thus 
have a tendency to become pioneers of the next genera- 
tion, such that they, too, enjoy abundant resources, just 
like their ancestors. Therefore, pioneer genes have a good 
chance to be carried along with the wave front and attain 
high frequencies, as if they "surf" on the wave. Thus, the 
descendents of an individual sampled from the tip of the 
wave have a finite probability to take over the wave front. 
In this case of "successful surfing" , further colonization 
will produce only descendents of the relevant pioneer be- 
cause the wave front has been "fixed" . The process of fix- 
ation at the front of a one-dimensional population wave 
is illustrated in Fig [I] 

The present study hinges on the question as to where 
lineages that reach fixation originate within the wave 
front. Clearly, the probability of successful surfing must 
increase with the proximity to the edge of the wave [l4| • 
On the other hand, more surfing attempts originate from 
the bulk of the wave where the population density is 
larger. We show that, due to this tradeoff, the origins of 
successful lineages have a bell-like distribution inside the 
wave front. Furthermore, this ancestral probability dis- 
tribution, together with the population-density profile of 
the population wave itself, is found to control the observ- 
able gradients in genetic diversity. The genetic pattern 
directly behind the moving colonization front turns out to 
mimic that of a small well-mixed (panmictic) population. 
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FIG. 1: Illustration of gene surfing by means of three consecutive snapshots of the genetic composition at the edge of an 
expanding population (in one spatial dimension), (a) A neutral red mutant arises at the wave front, (b) After some time, the 
genetic make-up at the wave front is drastically changed due to random number fluctuations and it is apparent that descendents 
of red will take over the wave front, (c) Fixation in the co-moving frame of descendents of red. Numbers in these sketches 
represent "inheritable" labels that are used in our simulations to trace back the spatial origin of individuals in the wave front. 
In this example, descendents of red are associated with position "9" in the co-moving frame. The dashed blue frame indicates 
the co-moving simulation box. 



The effective size N e of this population "bottleneck" is 
shown to be smaller than the typical number of individ- 
uals in the colonization front and very sensitive to the 
growth conditions in the very tip of the front. Coloniza- 
tion fronts in which individuals need to be accompanied 
by others in order to grow (Allee effect [HI) have a much 
larger effective population size than those in which indi- 
viduals grow even if they arc isolated from the rest. 

The outline of the paper is as follows. We first intro- 
duce a stochastic computer model that we use to gener- 
ate both pulled and pushed one-dimensional colonization 
waves. Tracer experiments within this model are then 
used to reveal the probability distribution of successful 
surfers and the decrease of genetic differentiation at the 
colonization front. Our succeeding theoretical treatment 
reveals to what extent both measures are related, and 
how they can be predicted for continuous models with 
quasi-stationary demography. After a comparison be- 
tween theory and simulations, wc discuss the significance 
of our results in the light of inferring past range expan- 
sions from spatially resolved genetic data. 



I. SIMULATIONS 

Most models for range expansions can be classified as 
describing pulled or pushed population fronts (l9l . [2p| . 
The distinction between the two cases corresponds to a 
difference in behavior. Suppose individuals need to be in 
proximity to other individuals in order to grow in num- 
ber (Allee effect [ill)- The presence of conspecifics can 
be beneficial due to numerous factors, such as predator 
dilution, antipredator vigilance, reduction of inbreeding 
and many others [2lj . Then, the individuals in the very 
tip of the front do not count so much, because the rate 
of reproduction decreases when the number density be- 
comes too small. Consequently, the front is pushed in 
the sense, that its time-evolution is determined by the 
behavior of an ensemble of individuals in the boundary 
region. On the other hand, a population in which an in- 
dividual reproduces, even if it is completely isolated from 
the rest, will be "spearheaded" by these front individu- 



als. These pulled fronts are responsive to small changes 
in the frontier and, therefore, are prone to large fluctua- 
tions 0. 

One might suspect that the genetic pattern left behind 
a population wave should reflect whether the colonization 
process is controlled by a small or large number of indi- 
viduals. Hence, we have set up a computer model that 
allows us to investigate the surfing dynamics for both 
classes of waves. 



A. Population dynamics 

The population is distributed on a one dimensional 
lattice, whose sites (demes) can carry at most iV individ- 
uals. The algorithm effectively treats individuals (•) and 
vacancies (o) as two types of particles, whose numbers 
must sum to N at each site. A computational time step 
consists of two parts: (i) a migration event, in which a 
randomly chosen particle exchanges place with a particle 
from a neighboring site. This step is independent of the 
involved particle types, (ii) A duplication attempt: Two 
particles are randomly chosen (with replacement) from 
the same lattice site. A duplicate of the first one replaces 
the second one (1 st — » 2 nd ) with probabilities based on 
their identities: proposed replacements o — > o, • — > • and 
• — ► o (growth) are realized with probability 1, whereas 
o — > • (death) is carried out only with probability 1 — s, 
depending on a growth parameter < s < 1. This asym- 
metry controls the effective local growth advantage of • 
over o. 

In terms of individuals and vacancies instead of parti- 
cles, we see that our model describes migration and local 
logistic growth of a population distributed over domes 
with carrying capacity N. Starting with a step-function 
initial condition, the simulation generates an expand- 
ing pulled population wave. The above algorithm rep- 
resents a discretized version [22j of the stochastic Fisher- 
Kolmogorov equation [23| with a Moran-type of breeding 
scheme [To| . To generate pushed waves as well, we extend 
our model by the following rule: In domes in which the 
number of individuals falls to N c or below, we set their 
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effective linear growth rate s to zero. This represents, for 
N c > 0, a simple version of the above mentioned Allee 
effect of a reduced growth rate when the population den- 
sity is too small. 

B. Tracer dynamics 

Tracer experiments within this computer model allow 
us to extract the genealogies of front individuals. After 
the population had enough time to relax into its prop- 
agating equilibrium state, all individuals are labeled ac- 
cording to their current position i G {1 . . .n} within the 
simulation box of length n, see Fig. QJi. These labels are 
henceforth inherited by the descendents, which thereby 
carry information about the spatial position of their an- 
cestors. The randomness in the reproduction and mi- 
gration processes (genetic drift) during the succeeding 
dynamics inevitably leads to a reduction in the diver- 
sity of labels present in the simulation box, see Fig. [TJd. 
Labels are lost due to either extinction or because they 
cannot keep up with the simulation box, which follows 
the propagating wave front Q. 

In our simulation, the gradual loss of diversity of labels 
at the wave front is measured by the quantity 

n 

ff(i)=J>(t)[l -!*(«)] > a) 

i=l 

which depends on the frequency Pi(t) of label i at time 
t after the wave has been labeled. H(t) represents the 
time-dependent probability that two individuals, ran- 
domly chosen from the bounded simulation box, carry 
different labels. Provided that mutations are negligible 
on the time-scale of the range expansion, we may think of 
our inheritable labels as being neutral genes at one par- 
ticular locus (alleles). We may thus identify H(t) with 
the probability that two alleles randomly chosen from the 
front region are different conditional on the well-mixed la- 
beling state at t = imposed by our simulation. Hence, 
we refer to H (t) as the time-dependent expected heterozy- 
gosity [To| at the wave front [45j . 

The perpetual loss of labels in our model without mu- 
tations eventually leads to the fixation of one label in 
the simulation box, see Fig. [TJ;. The value of this label 
indicates the origin within the co-moving frame of this 
successful "surfer" . It contributes one data point to the 
spatial distribution Pi of individuals whose descendents 
came to fixation. After fixation, the algorithm proceeds 
with the next labeling event. 

C. Results 

The parameters of our computer models are the deme 
size N, i.e. the maximal number of individuals per lattice 
site, the linear growth rate s per generation, and the crit- 
ical occupation number N c , below which the growth rate 



drops to zero (Alice effect). In our simulations, we set 
s = 0.1 throughout, and determine, for varying N and 
N c , the averages of the ancestral distribution function 
(Pi), the scaled occupation number (rii) /N, both being 
functions of the lattice site i in the co-moving frame, and 
the time-dependent probability of non-identity, (H(t)). 
Here, angle brackets indicate that the enclosed quanti- 
ties have been averaged in time, i.e. over many fixation 
events, and over multiple realizations of the same com- 
puter experiment [4^1 . 

Figure [2] illustrates the relation between the front pro- 
files (ni) /N and the ancestral distribution (Pi) in the co- 
moving frame. Whereas the wave profiles have the famil- 
iar sigmoidal shapes of reaction-diffusion waves 19, 20], 
the ancestral distribution functions are bell-curves with 
most of its support beyond the inflection point of the 
wave front. The fact that (Pi) has a maximum inside the 
wave front reflects a tradeoff, mentioned earlier, between 
a larger fixation probability in the tip of the wave ver- 
sus a larger number of surfing attempts originating from 
the bulk. Notice from Fig. [2^ that, for increasing deme 
size, the distribution becomes wider and shifts further 
into the tip of the wave, which is in contrast to the al- 
most TV— independent scaled wave profiles. Fig. [2)3 shows 
that the opposite effect is caused by increasing the cutoff 
value N c , which changes the type of the wave from pulled 
to pushed. 

Next, we measured the temporal decay of the heterozy- 
gosity H(t), defined in Eq. JT]). In Fig. [3J time-traces of 
H{t) are depicted for various parameters and show an ex- 
ponential decay after an initial transient. This allows us 
to characterize the strength of genetic drift at the wave 
front by a single number, the (asymptotic) exponential 
decay rate, — <9 t log (if (£)), which can be extracted from 
logarithmic plots of (H(t)). By analogy with well-mixed 
(panmictic) populations, in which the heterozygo sity de- 
cays exponentially with rate 2/N (Moran model[l(|), it 
is convenient to express the decay rate by 2/N e , in terms 
of an effective population size N e . The theoretical part 
below will further clarify to what extent the genetic di- 
versity at the wave front mimics that of a population 
"bottleneck" of constant size N e . 

Figure [4] depicts N e as a function of the deme size N 
on a double logarithmic scale for N c = and N c = 10. 
Naively, one might expect N e to be, roughly, the charac- 
teristic number of individuals in the width of the wave 
front, since these individuals contribute (by growing) to 
the advance of the wave. Thus, a linear relationship be- 
tween deme and effective population size would not be 
surprising. In contrast, we find that N e increases much 
slower than linearly with increasing deme size. Further- 
more, the effective population size turns out to be very 
sensitive to the presence of an Allee effect (N c > 0), 
which has the ability strongly increase the effective pop- 
ulation size. This point is illustrated, in particular, by 
the inset of Fig. Q] which depicts the effective population 
size N e in a simulation of fixed deme size (N = 1000) 
and varying strength of the Alice effect (10 < N c < 500). 
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(a) Pulled waves (N c = 0) (6) Pushed waves (N c > 0) 




i, lattice site i, lattice site 



FIG. 2: Measured distributions {Pi) (bell-curves) of "successful surfers" together with the normalized occupation numbers 
(rii) /N (sigmoidal curves; scaled along the vertical axis to fit the figure) as a function of the site number i in the co-moving 
frame; (a) for pulled waves (N c = 0) with varying deme sizes N; (&) for various pushed waves (N c > 0) with deme size TV = 1000 
compared to the corresponding pulled wave (dashed blue lines), which is also present in (a). 



(a) Pulled waves (N c = 0) 
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FIG. 3: The decay of genetic diversity, (H(t)), with time (in units of generations) on a log-linear scale for varying deme sizes; 
(a) for pulled waves (N c = 0); (b) for pushed waves with N c = 10. Both cases show, that an asymptotic exponential decay 
of (H(t)) is reached after an initial transient where the decay is weak. The duration of this transient is dependent on the 
size of the simulation box: The larger the simulation box, the larger the time until the exponential decay is approached. The 
asymptotic exponential decay rate, however, has been checked to approach a constant for a sufficiently large box size. This 
exponential decay rate is therefore well-defined and can be used to characterize the decrease of genetic diversity at the wave 
front. By analogy with panmictic populations, in which the heterozygosity decays exponentially with rate 2/N (Moran model), 
it is convenient to express the decay rate as 2/N e , i.e., in terms of an effective population size N e , which is noted in the legends, 
and plotted in Fig. [4] 



Qualitatively this phenomenon may be explained with 
the pushed nature of these waves. An Allee effect shifts 
the distribution P, of successful surfers away from the 
tip towards the wake of the wave (see Fig. [5b) and hence 
increases the gene pool from which the next generation 
of pioneers is sampled. This argument indicates a close 
relation between the N e and Pi, which also emerges ex- 
plicitely in the theoretical analysis below. 



II. THEORY 



The following employs a continuous reaction diffusion 
approach to establish a theoretical basis for the relation 
between the neutral genetic diversity and the popula- 
tion dynamics in non-equilibrium situations like range 
expansions, ft will help us to reconcile the somewhat sur- 
prising response of our simulations to parameter changes 
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FIG. 4: The measured effective population size N e as a 
function of deme size N on a log-log scale for pulled waves 
(N c = 0, asterisks) and pushed waves (N c = 10, crosses). The 
dashed and dotted lines have slope .30 and .42, respectively, 
i.e., significantly smaller than 1. Triangles represent the ef- 
fective population sizes as inferred from the strong-migration 
approximation, Eq. ©, using the measured (P)-distribution 
and population profiles. The inset shows the behavior of N e 
for varying cutoff-value N c and fixed deme size N = 1000, 
again, on a log- log scale. 



(deme size and Allee effect). Note from Fig. [2] that the 
changes in the ancestral distribution are dramatic, while 
the changes in the population profile itself arc quite mod- 
est. Results obtained from our approximation scheme are 
tested by direct comparison of simulations and theory. 



A. Gene surfing 

In our simulations, as well as in many other models 
of range expansions, a propagating population wave re- 
sults from the combination of random short-range migra- 
tion and logistic local growth. In the continuum limit, a 
general coarse-grained continuum description of such a 
reaction-diffusion system of a single species is given by 

d t c(x, t) = Dd 2 x c(x, t) + vd x c(x, t) + K(x, t) (2) 

formulated in the frame co- moving with velocity v, where 
c(x, t) represents the density of individuals at location x 
at time t and D is a diffusivity. The first two terms 
on the right hand side represent the conservative part of 
the population dynamics, for which we make the usual 



diffusion assumptions 2J|. The reaction term K(x,t) 



accounts for both deterministic and stochastic fluctu- 
ations in the number of individuals due to birth and 
death processes, and typically involves non-linearities 
such as a logistic interaction between individuals as well 
as noise caused by number fluctuations. For instance, our 
computer model with N c = maps, in the continuum 
limit [13], to the stochastic Fisher equation, for which 
K(x, t) = sc(coo — c) + ey/ c(coa — c)r), where r/(x, t) is a 



Gaussian white noise process in space and time, Coo oc N 
is the carrying capacity and e oc s/N sets the strength 
of the noise. We would like to stress, however, that the 
following analysis does not rely on a particular form of 
K. Therefore, we leave the reaction term unspecified. 

As in our tracer experiments, let us assume that inheri- 
table labels, representative of neutral genes, are attached 
to individuals within the population and ask: Given 
Eq. (J2j is a proper description of the population dy- 
namics, to what extent is the dynamics of these labels 
determined? To answer this question, it is convenient to 
adopt a retrospective view on the tracer dynamics. Imag- 
ine following the ancestral line of a single label located 
at x backwards in time to explore which spatial route its 
ancestors took. This backward-dynamics of a single line 
of descent will show drift and diffusion only] any reac- 
tion is absent because among all the individuals living 
at some earlier time there must be exactly one ancestor 
from which the chosen label has descended from. We 
may thus describe the ancestral process of a single lin- 
eage by the probability density G(£,r\x,t) that a label 
presently, at time t and located at x, has descended from 
an ancestor that lived at £ at the earlier time r. In this 
context, it is natural to choose the time as increasing to- 
wards the past, t > t, and to consider (£, r) and (x, t) as 
final and initial state of the ancestral trajectory, respec- 
tively. With this convention, the distribution G satisfies 
the initial condition G(£, t\x, t) = S(x — £), where S(x) is 
the Dirac delta function, and is normalized with respect 
to£, fG(£,T\x,t)d£ = l. 

Since t\x, t) as function of £ and r is a probabil- 
ity distribution function generated by a diffusion process 
that is continuous in space and time, we expect its dy- 
namics to be described by a generalized diffusion equa- 
tion (Fokker-Planck equation [HI)- Indeed, in the Ap- 
pendix El we show that, G(x,t\£,,r) obeys 

d T G(£,T\x,t) = -d e J(£,r\x,t) (3) 
J(£,r\x,t) = -DdtG + {v + 2Ddt\n[c(S,T)]}G , 

where all derivatives are taken with respect to the an- 
cestral coordinates (C,t). The drift term in Eq. ([3]) has 
two antagonistic parts. The first term, v, tends to push 
the lineage into the tip of the wave, and is simply a con- 
sequence of the moving frame of reference. The second 
term proportional to twice the gradient of the logarithm 
of the density is somewhat unusual. It accounts for the 
purely "entropical" fact that, since there is a forward- 
time flux of individuals diffusing from regions of high 
density to regions of low density, an ancestral line tends 
to drift into the wake of the wave where the density is 
higher. [47 

Our computer experiments measure the spatial distri- 
bution P of the individuals whose descendents came to 
fixation. This information is encoded in the long-time 
behavior of G 



P(f,r)= lim G(S,r\x,t) 



(4) 
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because it represents the probability that, in the far fu- 
ture (t — ► — oo, in our notation) when the population is 
fixed, an individual of the extant lineage has descended 
from an individual who lived at location £ of the co- 
moving frame at time r. Equation ((4]) has to be indepen- 
dent of x if fixation occurs: For t — > — oo, all individuals 
irrespective of their position x must have descended from 
the same ancestor and, thus, from the same location at 
the earlier time r. 

In principle, it is thus possible to relate the ances- 
tral distribution to the population dynamics by solving 
Eq. in the long-time limit. Unfortunately, this task is 
usually difficult to achieve analytically because the num- 
ber fluctuations in the density c(£, r) of the total popu- 
lation add noise to the drift term in Eq. ([3]). As is cus- 
tomary in many spatially explicit models of population 
genetics, let us suppose, however, that rules of "strict 
density regulation" 26| are imposed in order to guaran- 
tee a stationary demography, so that in the co-moving 
frame. 



c(x,t) « c st (x) . 



(5) 



Even though real systems and our discrete particle sim- 
ulations exhibit density fluctuations even in equilibrium, 
we take Eq. ((5J> as a first approximation in cases where 
the total number of particles is large enough, such that 
the relative magnitude of the density fluctuations is small 
(law of large numbers). We will call assumption Eq. ((Sj) 
the "deterministic approximation" as it neglects stochas- 
tic fluctuations in the total population density. 

With Eq. ((3J), all parameters in the Fokker-Planck 
equation, Eq. ||3J), of the ancestral distribution G(x, £|£, r) 
are time-independent and its analysis considerably sim- 
plifies: If a unique stationary solution P s t(£) = 
limt^oo t\x, t) exists, it can be written explicitely 
in terms of the stationary density profile c st (£), 



P st (£)cx Cs 2 t (Oexp«/D) 



(6) 



where a pre-factor is required to satisfy the normalization 
condition of P st , J P Bt (£, r\x, t) d£ = 

As shown below, the analytical expression Eq. (|6|) de- 
scribes at least qualitatively the bell-like shapes found for 
the ancestral distribution function (P,) in our stochastic 
simulations. The exponential factor biases the fixation 
probability [43] towards the tip of the wave (£ > 0) and 
competes with the pre-factor controlled by the decaying 
density of individuals in the tip of the wave. 

It is noteworthy that Eq. j6|) not only applies to range 
expansions, but can be evaluated for any deterministic 
population dynamics, such as deterministic models of 
evolution [U HH (where however rare events might be 
crucial as found in Ref. [2!|) and to source-sink popu- 
lations [13, [HJ, a simple example of which is given in 
the Appendix [B] If the spatial domain is unbounded, 
Eq. yields finite results as long as c s t(£) decays faster 
than exp[— v£/(2D)] as £ — > oo. This condition formally 
distinguishes the two classes of waves earlier denoted by 



pulled and pushed. Within the mean- field description of 
such waves, the right hand site of E q. (|6l) is normalizablc 
only in the case of pushed waves [32j. The density of 
pulled waves, however, decays as exp[— v£/(2D)] in the 
foot of the wave (£ — > oo) as follows from a linearized 
mean field treatment [33| . The prime example of pulled 
waves, the mean-field Fisher wave, does not therefore al- 
low for successful surfing, P st = 0, as is explicitely shown 
in Appendix [C] This is in marked contrast to our sim- 
ulations of stochastic Fisher waves (Fig. [2k)- There, we 
found finite bell-like ancestral distributions up to deme 
sizes on the order of 10 5 . This striking discrepancy in- 
dicates that the classical Fisher equation is a poor ap- 
proximation for the case of finite deme sizes (even if they 
are large). An improved deterministic equation with a 
modified reaction term has been proposed (34[, which is 
able to reproduce the leading reduction in wave velocity 
due to the discreteness. A remarkable property of Eq. ([6]) 
is that it should be valid, irrespective of the actual form 
of the reaction term, if the demography is determinis- 
tic. By comparing Eq. ([6]) to simulations, it is possible 
to test the deterministic character of a population wave, 
i.e., whether or not a deterministic reaction-diffusion de- 
scription might be appropriate. 

For our simulations, such a test is given in Fig. [5[ 
where we superimpose measured ancestral distribution 
functions (Pi) with those predicted by Eq. © based on 
the measured wave velocity v and the occupation num- 
bers (rii) (the discrete analog of the population density 
c s t(£)). It is seen that systematic deviations occur in 
the pulled case (N c = 0), where the predicted distribu- 
tion seems to be somewhat displaced towards the tip of 
the wave. The agreement of theory and simulation is 
much better for N c — 10 and further improves when N c 
is increased. Altogether, our deterministic approxima- 
tion Eq. ([6]) seems to apply best to pushed waves with 
a strong Allee effect (JV C 3> 1), whereas significant de- 
viations to Eq. ^ occur for pulled waves. An alterna- 
tive test of theory and simulation, presented in Fig. [6l 
supports this conclusion and furthermore shows that the 
deterministic approximation applied to pulled waves im- 
proves slowly with increasing deme size. For reasonable 
system sizes, however, fluctuation effects in pulled waves 
are non-negligible [l9[ . 



B. Decrease of genetic diversity 

To measure how fast genetic diversity decreases at the 
wave front due to gene surfing, we also studied in our 
simulations H(t), defined in Eq. fl} as the probability 
that two randomly sampled individuals carry different 
labels at a time t after a labeling event. To what extent 
are the decrease of H(t) and the shape of P(£) related? 
As before, a retrospective view on the problem simplifies 
the theoretical analysis. Imagine following the lineages 
of two randomly sampled individuals backward in time. 
They will drift and diffuse separately for a certain time t c 
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FIG. 5: The data from Fig. [3b (thin lines) superimposed 
with the prediction Eq. (JSJl (thick lines) for the ancestral dis- 
tribution function {Pi) based on the measured wave profile 
(rij) and velocity v. As explained in the text, the apparent 
systematic deviations in the pulled case (N c = 0) are caused 
by fluctuations in the tip of the wave. 



until they coalesce in the most recent common ancestor. 
If the last labeling event occurred at an earlier time t > 
t c (reversed time-direction) then both individuals must 
carry identical labels. If, on the other hand, t < t c then 
these individuals will have different labels unless their 
different ancestors happen to be in the same deme at the 
labeling time. Up to a small error of the order of the 
inverse size of the simulation box, we may thus identify 
the probability H(t) of two individuals carrying different 
labels with the probability that their coalescence time t c 
is larger than t. 

In principle, the coalescence time distribution of two 
lineages can be explored by studying the simultaneous 
backward-diffusion process of two lineages conditional 
on having not coalesced before [35|. This process is de- 
scribed by a generalization of Eq. ([2]) augmented by a 
well-known sink term [35| that accounts for the probabil- 
ity of coalescence when lineages meet. 

Here, we describe a (more tractable) approximation 
that estimates the behavior of H (t) from the distribution 
P(£), analyzed in the previous section. It is based on the 
assumption that the coalescence rate of two lineages is 
so small that each lineage has enough time to equilibrate 
its spatial distribution before coalescence occurs. Under 
this quasi-static approximation, the behavior of H(t) is 
described by [3f| 



-2H 



c(€,t) 



(7) 



when time is measured in units of generation times. The 
justification of Eq. ([7]) is as follows: The coalescence rate, 
—d T H, at time r in the past is given by the probability 
that the two lineages have not coalesced earlier, H(t), 
times the rate at which two separate lineages coalesce at 
time t. The latter is locally proportional to the product 



of the probabilities that the two lineages meet at the 
same place, oc P 2 (£,t), and that they meet in the same 
individual, cx c _1 , given they are at the same place. Less 
obvious, unfortunately, is the numerical pre- factor "2" 
on the right-hand side, which is specific to the employed 
breeding scheme (Moran model [lfj). 

Eq. ((7|) yields the correct coalescence time distribu- 
tion in the so-called strong migration limit [36l . \3l\ of 
large population densities, c — > oo, while the diffusiv- 
ity D and the spatial extension of the habitat are held 
fixed. In our case, it serves as a simple approximation 
that tends to overestimate coalescence rates, because it 
neglects spatial anti-correlations between non-coalescing 
lineages: Lineages that have avoided coalescence will usu- 
ally be found further apart than described by the product 
of (one-point) distribution functions in Eq. ([7p. Thus, 
their rate of coalescence will, typically, be smaller than 
in Eq. (0. 

Equation (J7J predicts exponential decay, H(t) ~ 
exp[— 2(t — t)/N e ], in the deterministic approximation, 
Eq. (f5l) , with a rate depending on a constant N c given 



n: 



j%(0 

Cst(0 



(8) 



In fact, a generalization of this argument to the coa- 
lescence process of a sample of n lineages shows that 
the standard coalescent [38( is obtained in the strong- 
migration limit with the parameter N e interpreted as the 
effective population size 36]. In other words, the coales- 
cence process in the strong-migration limit is identical, 
in every respect, to the coalescence of a well-mixed pop- 
ulation of fixed size N e . 

The strong-migration approximation may be tested by 
comparing the effective population sizes measured in our 
simulations with the ones predicted by Eq. © based on 
the measured ancestral distribution (Pj) and the number 
density profile {nf). These inferred values are plotted in 
Fig. [4] as (red) triangles. For both pushed and pulled 
waves, the agreement between inferred and measured ef- 
fective population sizes becomes excellent for large deme 
sizes, ./V > 100. For lower values of TV the strong migra- 
tion assumption overestimates genetic drift, presumably 
due to the neglect of correlations as mentioned above. 



III. DISCUSSION 

We have studied the impact of a range expansion on 
the genetic diversity of a population by means of simu- 
lations and analytical techniques. The one dimensional 
case treated in this article applies to populations follow- 
ing a (possibly curved) line, like a migration route, coast 
line, river or railway track. We have further simplified our 
analysis by neglecting habitat boundaries, which is ap- 
propriate for describing the colonization period, as long 
as the wave front is sufficiently far away from the bound- 
aries. 



8 



(a) Pulled waves (N c = 0) (&) fWia* waves (N c > 0) 




j, lattice site i, lattice site 



FIG. 6: The quantity ln((P») / (m) )/v as a function of lattice site i in the co-moving frame, which should have slope 1 according 
to the deterministic approximation, Eq. ((6]). (a) Results for pulled waves (N c = 0) with various deme sizes N; the values for the 
"initial slope" and "final slope" have been obtained by fitting a straight line to the lower and upper half of each shown curve, 
respectively, (b) Results for various pushed waves (N c > 0) with fixed deme size, N = 1000. (The domain of each curve is 
restricted to a region, in which the bell-like distributions (P;) has enough support to sample a sufficient amount of data points. 
This region roughly covers 98% of all successful surfing events.) 



Our findings suggest the following general scenario. 
Suppose, an initially well-mixed population increases its 
range from a smaller to a larger habitat, and that mu- 
tations may be neglected on the time scale of the range 
expansion. Our simulations show that the heterozygosity 
at the moving colonization front decays, due to genetic 
drift, exponentially in time with a rate 2/N e , depend- 
ing on an effective population size N e . Upon combining 
this rate with the velocity v (per generation) of the colo- 
nization front, we obtain a length A = vN e /2 that char- 
acterizes the pattern of genetic diversity generated by 
the colonization process: As the wave front moves along 
it leaves behind saturated denies with heterozygosities 
given by the value of the front heterozygosity as the wave 
front passes through. The transient colonization process 
therefore engraves a spatially decreasing profile of het- 
erozygosity into the newly founded habitat. This profile 
decays exponentially in space on the characteristic length 
A and serves as an initial condition for the succeeding 
period of demographic equilibration, which may be de- 
scribed by traditional models of population genetics [l(| • 

Of critical importance for the interpretations of gene 
frequency clines in natural populations is the question 
as to which parameters chiefly control A. Our computer 
simulations have revealed, that the effective population 
size, and thus A, only grows sub-linearly with increas- 
ing deme size, in contrast to the naive expectation that 
the effective population size should roughly be given by 
the characteristic number of individuals contained in the 
wave front (deme size times the width of the wave). On 
the other hand, we found that the population "bottle- 
neck" at the wave front was significantly widened when 
we implemented an Allee effect [l8T | into our model, by 



which growth rates for small population densities arc de- 
creased. As a consequence, the region of major growth 
shifted away from the frontier into the bulk of the wave 
and, thus, the effective size of the gene pool for the fur- 
ther colonization was increased. 

Finally, we have developed a theoretical framework to 
study the backward-dynamics of neutral genetic markers 
for a given non-equilibrium population dynamics, which 
is summarized by the Fokkcr-Planck equation Q and its 
generalization in Appendix lAl In stationary populations, 
it leads to a simple expression, Eq. ^ , for the long-time 
probability distribution of the common ancestors, which 
can be used, in the strong-migration limit to determine 
the effective population size, via Eq. |(5J). Comparison 
with our stochastic simulations reveals that the simple 
deterministic results are good approximations to stochas- 
tic simulations in the case of pushed waves (strong Allee 
effect), but that significant deviations occur for pulled 
waves due to the fluctuations in the frontier of the pop- 
ulation wave. For sufficiently large deme sizes, the effec- 
tive population size could, in both cases, be inferred with 
remarkable accuracy from Eq. ([5]). 
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APPENDIX A: BASIC FORMALISM 

In this appendix, we describe the derivation of the 
Fokkcr-Planck equation, Eq. ([3]), and its generalization 
to heterogeneous migration and higher dimensions. 

Central to our analysis is the assumption that the total 
population consists of statistically identical entities, such 
that two different individuals at a given time and location 
behave in the same way. In particular, we assume that 
migration as well as reproduction of an individual are 
exchangeable [391 ] random processes, i.e., independent of 
any label that might be assigned to it. 

Let us indeed imagine a subpopulation labeled by a 
neutral marker, and ask: What is the dynamics of the la- 
beled individuals for a given dynamics of the total popu- 
lation, as described by Eq. (fj))? It is clear that the density 
c*(x, t) of labels obeys a reaction-diffusion equation with 
coefficients D*, v* and K* that are closely related to the 
one for total population by means of the exchangeability 
assumption. Firstly, labeled and unlabeled individuals 
should migrate statistically in the same way, which is 
measured by diffusion and drift coefficients, i.e., we have 
D* = D and v* = v. Secondly, the labeled subpopu- 
lation must carry the number fluctuations of the total 
population, encoded in K, in proportion to its reduced 
size: K* = Kc*/c. However, these statements are true 
only on average: The discreteness of the particle num- 
bers lead to fluctuations, for instance, in the quantity 
K* — (c* /c)K. To illustrate this point, imagine that at a 
given time, the term K dictates that an individual dies in 
some small spatio-temporal region, then this individual 
has to be sampled from the labeled subpopulation with 
probability c*/c. The fluctuations of c* due to this sam- 
pling procedure represents random genetic drift and must 
have zero mean according to the exchangeability assump- 
tion. Similar fluctuations affect the migration currents 
of the labeled subpopulation. Thus, only upon averag- 
ing over this source of stochasticity, we may formulate a 
reaction-diffusion equation of the form 

d t c* = Ddlc* + vd x c* + K C — ( Al ) 

for the average density c*(x, t) of labeled individuals. In 
Eq. (|AI[) , only the effect of the genetic drift of labeled in- 
dividuals within the total population has been averaged 
out, the number fluctuations affecting the total popula- 
tion density are retained through the fluctuating reaction 
term K. In other words, Eq. (|A1[) describes the behavior 
of labeled individuals, if we average over many realiza- 
tions conditional on a given fixed evolution of the total 
population, described by Eq. d2|). Note that, this aver- 
aging "works" because Eq. (|A1[) is linear in c*, such that 
a noise term with zero mean added to the right-hand- 
side to account for the genetic drift can be averaged out 
without generating higher moments. 
Next, we use 

K(x, t) = d t c(x, t) — Dd%.c(x, t) — vd x c(x, t) 



from Eq. ([2]) to substitute the reaction term K in 
Eq. (jAip . After rewriting the average density c*(x, t) = 
p{x, t)c(x, t) of labeled individuals in terms of their aver- 
age frequency (or ratio) p, we obtain 

d t p(x, t) = DB%p + {v + 2Dd x ln[c(x, r)]} d x p . (A2) 

Notice that, p =const. is a (steady state) solution of 
Eq. (|A2[) . Relations formally equivalent to Eq. (|A2|) have 
been formulated in Refs. [4CJ,[4JJ in different contexts, in 
which, unlike the present case, random genetic drift could 
not be averaged out due to non-linearities, but had to be 
disregarded, instead. As a transport equation for deter- 
ministic gene frequencies, an equation similar to Eq. (|A2[) 
was also recently obtained in Ref. [25[. 

For a given realization of the time-evolution of the to- 
tal population, the quantity p(x, t) can be interpreted 
as the probability that an individual sampled from (x,t) 
is labeled, i.e., that it has descended from the initially 
labeled population. If the above tracer experiment for 
a given dynamics of c(x, t) is repeated multiple times, 
p(x, t) represents the histogram of the number of times a 
individual sampled from (x, t) is labeled. 

By choosing proper initial conditions, this allows us 
to study "where individuals come from": Suppose that, 
at time r, the labeled population contains all individu- 
als within a small interval around the position £. The 
solution of Eq. (|A2[) for later times will then tell us the 
probability that a individual at (x, t) has descended from 
an ancestor sampled from that narrow region around £ 
at time r. 

Hence, the probability density G(£, t\x, i), introduced 
above Eq. ([3]), that an individual at (x, t) has an ancestor 
who lived at (£, r), is just the solution of Eq. (|A2[) for the 
initial condition 

G(£,0\x,0) = 6(x - $ , 

where S(x) is the Dirac delta function. It is straightfor- 
ward [24[ to show that this Green's function of Eq. (|A2[) 
also obeys the Fokker-Planck equation ([3]), in which time 
is measured in the backward direction; Eq. (IA2I ) is usually 
called the (Kolmogorov-) backward equation [24[ associ- 
ated with the Fokker-Planck equation ©. 

The content of the Fokker-Planck equation ([3]) may be 
further illustrated by a physical analog. The dynamics of 
a single ancestral line backward in time conditional on a 
particular demographic history, as described by Eq. l[2]). 
is equivalent to the Brownian motion of a particle (with 
fcfiT = 1) in a potential E/(£, t), whose negative gradient 
is given by the drift coefficient v + 2Dd^ ln(c) in Eq. ([3]). 
The form of this time-dependent potential, 

U(^r) = -v^-2Dhx[c^,T)} , (A3) 

suggests an interpretation in terms of a fluctuating free 
energy in which the first and second term are the en- 
ergetic and entropic contribution, respectively. If the 
spatial domain is unbounded, the long-time distribution 
function of the fluctuating particle will be non-trivial only 
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if the potential Eq. (j A3|) has the form of a well, which 
is able to "trap" the fluctuating particle on long times. 
Otherwise, for instance if the potential is half-open as 
in the case of a Fisher wave, the ancestral probability 
distribution will decay to zero on long times. 

So far, our analysis was restricted to a one-dimensional 
reaction diffusion model, in which drift and diffusion are 
space-independent. The generalization of Eq. @ to het- 
erogeneous migration and higher dimensions reads 



dc 







(A4) 



The conservative dynamics is now described by a ma- 
trix Dijix, t) of diffusivities and a velocity vector Vi(x, i), 
which may both depend on space and time. Again, 
the reaction term K(x, t) accounts for deterministic and 
stochastic fluctuations in the number of individuals due 
to birth and death processes. 

Repeating the above arguments for the dynamics of 
neutral labels under the more general population dynam- 
ics Eq. (|A4[) yields a multi-dimensional Fokker-Planck 
equation for the probability density G(£, t\x, t) that an 
individual at (5? , t) has descended from an ancestor who 
lived at (£,t), which is given by 



d T G(£,T\x,t) 

Ji(£,T\x, t) 



d 

-|-(A^) + 



(A5) 



3D 



at 



, ± ^ 2 d{D lj c) 
c d$) 



G . 



Here, summation over identical indices is implied and 
time again increases in the backward direction. The nat- 
ural requirement that there is no probability flux J out 
of the region S of non-vanishing population density leads 
to reflecting boundary condition, J = 0, on the bound- 
ary dS of S. Note that our stochastic description of the 
backward dynamics of a single lineage, Eq. (|A5|) . is fully 
determined by the demographic history c(x,t). A knowl- 
edge of the actual form of the reaction term K(x, t) in 
the reaction-diffusion equation (|A4[) is not necessary. 



APPENDIX B: SOURCE-SINK POPULATIONS 

For purely conservative populations [42| of neutral in- 
dividuals, subject only to diffusion and drift, it is well- 
known that the fixation probability is the same for all 
individuals, u =const.= 1/N, and that the effective pop- 
ulation size equals the total population, N e = N. How- 
ever, when reaction terms are important, individuals be- 
come privileged or handicapped depending on where they 
linger. As a telling example, let us consider the case of 
an "oasis" [3l| (or source [43|) with a carrying capacity 
c in equilibrium with a "desert" (a sink) of smaller car- 
rying capacity Cd- Sufficiently far away [50( from the con- 
tact zone, the population densities will be saturated at 



their respective carrying capacities. According to Eq. ([6]), 
the ancestral distribution function of a stationary non- 
moving population will be locally proportional to the 
square of the density, -P s t(£) oc c^ t (£). As noted in foot- 
note 5 of the main text, the fixation probability of a 
mutation occurring in a single individual is given by the 
ratio of ancestral distribution function and population 
density, «(£) oc P s t(£)/c s t(£). Thus, the probability that 
a neutral mutation fixes will be larger if it arises deep in 
the oasis, it , than if it arises in the desert, Ud- The ratio 
of both fixation probabilities is given by the ratio of the 
respective densities, 



Uo 
Ud 



Co 
Cd 



(Bl) 



Apart from its simplicity, the relation Eq. (|B1[) is remark- 
able because it is independent of the diffusion constants 
and the details of the particular logistic interaction be- 
tween individuals. 

If the combined system of oasis and desert is closed, 
the effective population size, Eq. ([8]), evaluates to 



(c 2 L d + c 2 L f 
c\L d + c%L a 



< CdL d 



(B2) 



to leading order in the linear sizes L a and Ld of oasis and 
desert. As mentioned earlier, our reasoning regarding 
coalescence times only applies to the strong-migration 
limit, in which the fixation time ~ N e T, where T is 
the generation time, is much smaller than the longest 
relaxation time of the Fokker-Planck equation. For the 
present case, the latter may be estimated by the time 
needed for lineage to cross the habitat, ~ L 2 /D. 



APPENDIX C: SURFING ON A FISHER WAVE 

Here, we apply our theory of gene surfing to the Fisher 
equation (23j . 



d t c(x, t) = Dd 2 x c(x, t) + s 



c(x, t) 



c(x,t), (CI) 



which is the prime example of a pulled front. This equa- 
tion was originally proposed as a mean-field model for 
the spread of a dominant gene with selective fitness ad- 
vantage s through a population with constant density Coo 
(carrying capacity). Equation (|C1|) has also been useful 
as a description of an expanding population, for which s 
is the difference between the linear birth and death rates, 
and the term — sc 2 /coo represents some self-limiting pro- 
cess, roughly proportional to the number of pairs of in- 
dividuals at position x. 

There are two spatially homogeneous fixed points: an 
unstable fixed point at c(x) = 0, in which there is no 
population at all, and a stable fixed point at c(x) = c 001 
where the population saturates to the carrying capacity 
Coo of the environment. 
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Non-negative initial configurations evolve smoothly to- 
ward the stable fixed point; analysis of the time devel- 
opment of spatial fluctuations in this model reveals that 
equilibrium can be reached via traveling soliton-like so- 
lutions Cgt (x — vt) , referred to as Fisher waves [23| , which 
represent steady state solutions of Eq. (|Clj) . In the wave 
front, where the population density is much smaller than 
the carrying capacity, the non-linear logistic term oc c 2 in 
Eq. (|C1[) may be neglected. The steady state solution of 
the remaining linear equation is exponentially decaying, 

c st ~ exp(-x/\) , (C2) 

for x — > oo, where decay length A and velocity v are 
related by 

D v 

° = ^-x + s - (C3) 

The population density to be nonnegative requires real 
solutions A > of Eq. (|C3|) , which do not exist unless 

v > 2VDs . (C4) 

It can be shown that the solution corresponding to the 
lowest velocity v = 2\J Ds and decay length A = -J D / ' s is 
approached for any initial conditions with compact sup- 
port, and, thus, is the solution most relevant for biologi- 
cal applications [33j]. 

Now, when we evaluate the surfing probability, Eq. ^ 
for this standard model of a spreading wave, we find zero - 
a somewhat surprising result in light of the finite bell-like 
curves obtained from our simulations of stochastic Fisher 
waves (Fig. 2a). The exponential decrease of the popu- 
lation density at the wave front, Eqs. (|C2[|C4p . is simply 
not fast enough to render the function c 2 t (£) exp(u£/D) 
normalizablc - not even for the lowest velocity for which 
it approaches a positive constant as £ — * oo. Hence, a 
non-zero stationary distribution function of the common 
ancestor does not exist, even though the total population 
is in a steady state. 



A closer look to the dynamics of a lineage, as de- 
scribed by the Fokker-Planck equation reveals how 
the distribution of the common ancestor decays to zero 
with time. As we evolve the probability distribution 
G(£,t\x, t = 0) that a lineage diffuses from a location x 
to £ back through time, it spreads out, due to diffusion, 
and is subject to a drift of strength v + 2Dd^ logc s t- If a 
lineage starts out deep in the wake of the wave, x <C — 1, 
where logc s t — > 0, it experiences a drift pushing it to- 
wards the wave front. After a time of the order of |x| /v, 
the probability cloud of the single lineage reaches the 
front and, when the inflection point is passed, the drift 
has decreased appreciable because — c\logc st = 0(1). In 
the tip of the wave, the density profile is exponential, 
Eq. (|C2|) . such that the drift, v + 2Dd^ log(c st ), saturates 
at a non-negative value w = (v — 2D/X) > 0, which 
tends to move the lineage even further away from the 
bulk into the tip of the wave. For large times, the distri- 
bution G assumes the form of a bell-like curve of width 
~ s/Dt moving with a velocity w. For any fixed £ in 
the foot of the wave and w > 0, the distribution func- 
tion G decays exponentially to zero for times, when the 
probability "cloud" has passed £ . Only for w = 0, which 
corresponds to the lowest allowed velocity v = 2\J sD, 
drift is absent in the tip of the wave, such that the decay 
is much slower, G — > t -1 / 2 . In any case, however, G 
decays to zero for any location £, which means that, no 
matter which individual we choose, the fixation proba- 
bility is zero. Thus, "successful surfing" is not possible 
in the case of the deterministic Fisher wave, as was men- 
tioned already in the main text. We expect that, in a 
stochastic simulation of a finite number of particles, the 
time-dependent features discussed above for the mean- 
field Fisher equation are merely transient and only visible 
for times smaller than the longest relaxation time of the 
Fokker-Planck equation, Eq. ([3]). This conjecture can be 
supported by employing an approximation scheme due 
to Brunet-Derrida [34| to take into account leading order 
effects of discreteness. 
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